Projet: Contrôle optimal par recherche arborescente et algorithme UCT

Matthieu Guerquin-Kern1

Le 16 janvier 2007
Abstract: Disposant d'un modèle de la dynamique d'un processus, on peut envisager d'aborder la résolution du problème de contrôle par une recherche arborescente. L'arbre est alors construit à partir de diverses simulations Monte-Carlo. Dans ce travail, l'application visée est le contrôle du processus de digestion anaérobie de matière polluante organique par une colonie de bactéries. Nous disposons d'un modèle du processus, mais il est relativement complexe. Notre arbre de recherche arborescent se construit donc dans un espace d'états à grandes dimensions, ce qui rend rapidement incommensurable le nombre de calculs à effectuer. C'est pourquoi nous envisageons l'utilisation du récent algorithme, UCT, qui, résolvant le fameux dilemme exploration-exploitation, optimise le parcours de l'arbre dans un temps donné.

Contents

1  Introduction

Le cahier des charges de ce travail est de fournir une méthode de contrôle permettant d'optimiser une production dans le cadre d'un processus de digestion anaérobie. Disposant de modèles dynamiques de ce processus [1], on peut envisager de simuler tous les scénarii possibles à partir d'un état donné et de choisir l'action qui permet dans le futur de maximiser la production. Cet algorithme retourne alors une action qu'on souhaite proche de l'optimum si la profondeur d'exploration est suffisante.

Ainsi, à chaque fois que l'on cherche l'action optimale à l'instant t, on va construire, par diverses simulations, un arbre dont les noeuds sont les états et dont chaque action est représentée par une branche menant vers un nouvel état. L'action choisie sera alors la première action de la séquence qui mène au bout de l'arbre à la production maximale. Bien entendu, cet algorithme nécessite d'avoir un ensemble d'actions possibles de cardinal discret. Autrement dit, de quantifier les actions. Ainsi, l'on comprend que si on explore l'arbre suffisamment loin et avec un ensemble d'action suffisamment grand, on pourra très bien approcher l'action optimale.

Le problème qui se pose alors est celui de la construction de l'arbre. Dans un premier temps, pour valider la méthode, nous construirons l'arbre en quantifiant les actions possibles et en explorant la totalité des séquences d'actions dans A et ce jusqu'à une certaine profondeur. Mais alors, le nombre d'états à explorer devient rapidement gigantesque. Ce qui pose le problème de la faisabilité des calculs en temps réel. En effet, tous ces calculs seront limités par la fréquence à laquelle on se pose la question de l'action optimale. Toutefois, dans le cas de la digestion anaérobie, les constantes de temps peuvent être longues donc la fréquence d'échantillonnage relativement réduite.

Dans un second temps, on cherchera à parcourir l'arbre de façon optimisée, afin que, dans le temps imparti, on ait perdu le moins de temps de calcul possible dans des scénarii vains et que par contre on ait exploré au maximum les voies les plus prometteuses. Autrement dit, on souhaiterais passer utiliser au maximum le temps de calcul disponible dans l'exploitation au détriment de l'exploration. Mais, on sait aussi qu'on ne peut pas totalement négliger l'exploration. Ces deux nécessités antagonistes sont pourvues pour le problème du bandit à armes multiples par un algorithme appelé UCB [2] (pour Upper Confidence Bounds). Récemment, cette solution a été adaptée au cas de la recherche arborescente dans un algorithme appelé UCT [3] (pour Upper Confidence bounds in Tree). Nous chercherons donc à mettre à profit l'efficacité de cet algorithme pour notre application.

2  Description du problème

2.1  Le contexte

On envisage de traiter les eaux usées par un processus de digestion anaérobie. Ce procédé est basé sur un écosystème complexe d'espèces bactériennes anaérobies dégradant la matière organique. Par rapport aux techniques de traitement aérobies plus classiques, elle peut dégrader des substrats difficiles à des concentrations élevées, ne produit que très peu de boues, nécessite un minimum d'énergie et peut même produire de l'énergie en utilisant la combustion du méthane (co-génération) et rapidement devenir rentable. La contrepartie est que ce procédé est très délicat à contrôler. En effet, une perturbation, même petite, peut conduire à la déstabilisation du procédé du fait de l'accumulation de composés intermédiaires toxiques au sein du réacteur conduisant à l'élimination de la biomasse2.

2.2  Modèles du processus

2.2.1  Modèle simplifié

Nous utiliserons ce modèle simplifié au départ, afin de réduire la dimension de l'espace d'état (ici 2).
St+1−St
τ
=
at(Sin − St) − k1Xt
c1
c2St+c3/St
    (1)
Xt+1−Xt
τ
=
−α atXt + k2Xt
c1
c2St+c3/St
    (2)
Ot+1 =
k3Xt
c1
c2St+c3/St
    (3)
avec at le débit entrant en substrat, S la concentration en substrat, Sin la concentration en substrat entrant, X la population bactérienne, O la production de biogaz, τ le pas de temps, ci les constantes de réaction et ki, α des constantes.

Le terme at(Sin − St) traduit la dilution ou la concentration en substrat. Le terme − k1Xtc1/c2St+c3/St représente la consommation de substrat par les bactéries. Le terme k2Xtc1/c2St+c3/St figure la croissance de la population bactérienne. Le terme −α atXt traduit les pertes en bactéries du fait de la dilution.

Nous prendrons les valeurs suivantes: τ=1, Sin=0.5, α=0.1, c1=1, c2=3, c3=8, k1=1, k2=2 et k3=2.

2.2.2  Modèle plus élaboré

Ce modèle est celui décrit dans [1]. La dimension de l'espace d'état est 4.

S1,t+1−S1,t
τ
= at(S1,in − S1,t) − k1μ1(.)X1,t     (4)
X1,t+1−X1,t
τ
= −α atX1,t + μ1(.)X1,t     (5)
S2,t+1−S2,t
τ
= at(S2,in − S2,t) +k2μ1(.)X1,t − k3μ2(.)X2,t     (6)
X2,t+1−X2,t
τ
= −α atX2,t + μ2(.)X2,t     (7)
Ot+1 = k4μ2(.)X2,t     (8)
avec Si les concentrations en substrats 1 et 2, Si,in les concentrations en substrat entrant. Xi les populations bactériennes 1 et 2, O la production de biogaz, τ le pas de temps et ki, α des constantes.

Les μi(.) sont les facteurs de croissance pour les populations 1 et 2. On supposera qu'ils suivent les lois: μ1(.)= c1S1,t/c2+S1,t (cinétique de Monod) et μ2(.)= c3/c4S1,t+c5/S1,t (cinétique de Haldane).

Nous prendrons les valeurs suivantes: τ=1, S1,in=0.5, S2,in=0.1, α=0.1, c1=1, c2=8, c3=1, c4=3, c5=8, k1=1, k2=2, k3=2 et k4=2..

2.3  Processus de décision markovien

Notre problème peut se voir comme un problème de décision Markovien.

2.4  Formalisation

Si on considérait le processus comme étant stochastique, on pourrait représenter un embranchement de l'arbre comme en 1(a)5Subfigure 1(a)subfigure.1.1.
[cas stochastique]    [cas déterministe]

Figure 1: embranchement de l'arbre


Alors, la valeur Vx à attribuer à un noeud x de l'arbre serait le maximum des valeurs V(x,a) attribuées aux états-actions (x,a) lui étant liés: Vx=maxaV(x,a). Ces valeurs-ci étant elles-mêmes estimées en moyennant les récompenses r obtenues dans la transition du noeud état-action vers le noeud état suivant et en l'ajoutant à la valeur des états suivants x': V(x,a)=E(r+Vx'|a). On attribue la valeur 0 aux états ”feuilles” de l'arbre et par descente vers la racine, on attribue les valeurs à chaque état. Ainsi, à l'état racine, l'action choisie sera celle pour laquelle la valeur du couple état-action est la plus grande.

Dans le cas déterministe, les choses sont légèrement plus simples. On peut représenter un embranchement comme en ??. Il n'est alors plus besoin de faire intervenir explicitement les noeuds état-action. Nous cherchons la trajectoire dans l'arbre qui possède la plus grande somme de récompenses actualisée. L'idée est alors de ne plus calculer la valeur de chaque état mais directement la somme actualisée des récompenses obtenues sur le trajet depuis la racine jusqu'à chaque noeud. Cette somme peut se mettre à jour directement pendant la construction de l'arbre. Dès lors, si on veut estimer la valeur d'un noeud de l'arbre, il suffit de prendre la somme actualisée maximale sur les feuilles de l'arbre issues de ce noeud. Une fois l'arbre calculé, l'action à choisir sera la première action de la trajectoire partant de la racine et menant à la feuille de somme de récompenses actualisée maximale.

Si on se donne M actions possibles et que l'on explore entièrement l'arbre sur une profondeur p, on aura Σi=1pMi états à calculer et Mp feuilles à l'arbre.

3  Recherche arborescente complète: une approche naïve

Dans cette partie, nous testons le principe du contrôle par recherche arborescente. Nous cherchons également à caractériser l'influence des paramètres que sont la profondeur d'exploration p, le nombre d'actions possibles M ou encore le coefficient d'actualisation γ. Le modèle 2.2.14Modèle simplifiésubsubsection.2.2.1 du procédé est utilisé car cette méthode est coûteuse en calculs.

Les figures 26Recherche arborescente complètefigure.2 et 37Recherche arborescente complètefigure.3 présentent la construction des arbres.


Figure 2: Représentation d'un arbre dans l'espace d'état avec 2 actions possibles sur une profondeur de 5.




Figure 3: Représentation d'un arbre dans l'espace d'état avec 3 actions possibles sur une profondeur de 4.


On a représenté le point de départ en rouge et les feuilles sont représentées par un point vert.

Á partir du moment où la construction de l'arborescence depuis un état racine avec la donnée des actions possibles fonctionne, nous pouvons mettre en oeuvre l'algorithme de pilotage du procédé.

La figure ?? illustre ce contrôle par l'évolution de l'état sur 100 itérations avec 3 actions possibles et une profondeur d'exploration de 5. Le coefficient d'actualisation est choisi à γ=0.9.


Figure 4: Fonctionnement de l'algorithme sur 100 itérations


On remarque qu'au départ, l'action est nulle. La population bactérienne grandit et consomme l'excès de substrat. C'est le comportement du système lorsqu'il est laissé au repos (pas d'alimentation en substrat). Puis la population se régule et on a convergence de l'état vers un régime permanent.

3.1  Influence du coefficient d'actualisation

Étant donné le problème, si on veut accorder une importance aux états futurs, il faut prendre γ proche de 1 (par valeur inférieure, sinon on devient instable). On pourra mettre en relation cette grandeur à l'incertitude que l'on peut avoir sur la prédiction d'un état futur à un horizon temporel donné. La valeur donnée à γ doit également être considérée avec la profondeur d'exploration de l'arbre. En effet, si γ est faible devant 1, quelque soit la profondeur d'exploration, les récompenses auront rapidement avec la profondeur un poids très faible.

On garde la même profondeur d'exploration (5) et le même nombre d'action (3) pour pouvoir comparer les résultats avec le test précédent, en faisant varier le coefficient d'actualisation.

On constate que plus γ est faible, plus on converge lentement et régulièrement vers un équilibre asymptotique. Cependant, plus γ est faible, plus la valeur moyenne de la production en cet état est faible. Prendre γ faible n'est donc pas, comme prévu, une bonne solution.

Si on cherche à prendre γ le plus proche possible de 1 (par valeur inférieure), plus l'amplitude des oscillations autour du régime permanent sont grandes. Par contre la rapidité de convergence est toujours très bonne et la valeur moyenne de la production s'améliore un peu.

Si on prend le cas limite, γ=1, les oscillations restent contrôlées car en fait la profondeur d'exploration est faible. Si on prend γ=1.1, on observe bien une croissance des oscillations.

3.2  Influence du nombre d'actions

On s'attend à ce que le nombre d'action augmente le domaine d'exploration de l'algorithme et donc son efficacité. Mais, du fait de son lien avec le nombre de noeud dans l'arbre d'exploration (Σi=1pMi) et donc la rapidité de calcul de l'algorithme, ce degré de liberté est à utiliser avec grande parcimonie.

On constate, toujours dans les conditions de l'expérience de la figure ?? et avec un nombre d'actions possibles de 2, 3 puis 4, que l'état moyen en régime permanent ne change guère. Par contre, quand on augmente le nombre d'actions possibles, l'amplitude des oscillations du régime permanent diminue. Cela se comprend aisément si on suppose qu'il existe une valeur d'action pour laquelle, maintenue en permanence, l'état reste stationnaire. Alors, en régime permanent, l'action quantifiée que l'on va appliquer devra varier pour atteindre en moyenne cette valeur stationnaire. Plus fine est la quantification de l'action (plus grand est M le nombre d'actions possibles), plus faibles seront les variations de l'action quantifiée.

3.3  Influence de la profondeur de l'exploration

On retrouve le même comportement de l'algorithme avec la profondeur qu'avec le coefficient d'actualisation. C'est-à-dire que pour de faibles valeurs la convergence est plus lente et la valeur moyenne de production est plus faible mais le régime permanent est plus régulier. Ceci conforte l'analyse que l'on a pu faire.

En outre, dans le cas où la profondeur est seulement de 1, on a aucune vision du problème à long terme et la politique suivie consiste à ne jamais alimenter en substrat, ce qui, certes, permet d'augmenter à très court terme la population et donc la production mais l'épuisement du substrat inverse vite la tendance.

3.4  Critique

On a vu avec la figure ?? que le contrôle du processus par ce résultat peut donner des résultats correct. Mais, la figure 59Critiquefigure.5 montre bien que cette méthode de contrôle est ridicule (en tout cas pour cette profondeur d'exploration). En effet, dans les mêmes conditions d'expérimentation, une loi d'action constante donne de loin de meilleurs résultats.


Figure 5: Fonctionnement de l'algorithme avec une action constante.


D'ailleurs, remarquons que, le modèle 2.2.14Modèle simplifiésubsubsection.2.2.1 étant relativement simple, on pourrait prédire la valeur constante à donner à l'action pour maximiser la production. Dès lors, notre problématique perd beaucoup de son intérêt.

4  Construction de l'arbre par UCT

On a pu constater avec le travail fait précédemment que pour des raisons de temps de calcul, on est rapidement limité en dimension de l'espace d'état, en profondeur d'exploration et en nombre d'actions. Or, pour un bon contrôle, les deux derniers paramètres doivent être suffisamment grands. On cherche donc à utiliser une méthode de parcours d'arbre plus optimisée.

Le principe décrit dans cette partie est fortement inspiré de [3]. C'est celui que j'ai mis en oeuvre dans le code Matlab qui se trouve en Annexe.

Les deux principales caractéristiques voulues pour cet algorithme sont
(i)
Faible probabilité d'erreur si l'algorithme est stoppé prématurément,
(ii)
Convergence vers l'action optimale si assez de temps est fourni.

4.1  Motivations

On a vu dans la partie précédente que, dans l'arbre, les états visités pouvaient être très proches. D'où, dans un premier temps, l'idée de quantifier les états afin de ne pas perdre du temps de calcul dans des morceaux de trajectoires quasi-redondants.

Remarque 1: Le fait de quantifier les état prédits engendre une imprécision supplémentaire du contrôle. On veillera donc à avoir des pas de quantification suffisamment petits. En outre, on veillera à mettre la quantification en accord avec la quantification des actions (en ce sens, qu'une quantification trop grossière des états pourrait rendre superfétatoire la finesse de quantification des actions).

Remarque 2: On pourrait également songer à ne pas quantifier les états systématiquement mais simplement définir un critère de voisinage avec les autres états de l'arbre et éventuellement faire coller le nouvel état avec un de ceux existants si le critère de voisinage est remplit.

Les pas de quantification sont finalement choisis grâce aux dimensions caractéristiques des ”paquets” d'états mis en évidence dans les figures 26Recherche arborescente complètefigure.2 et 37Recherche arborescente complètefigure.3.

D'autre part, on espère pouvoir gagner en rapidité de parcours si, au lieu de visiter l'intégralité de l'arbre, on se concentre sur les états les plus prometteurs (exploitation) sans toutefois négliger complètement l'exploration.

4.2  Description de l'algorithme UCT

On va construire l'arbre en lançant un nombre fixé de trajectoires à partir de la racine.

Sur une trajectoire, pour chaque noeud, on se pose les questions suivantes: Cet algorithme fait apparaître quatre degrés de liberté:
  1. le nombre de trajectoires lancées , sur lequel on pourra jouer pour changer la taille de l'arbre.
  2. la profondeur maximale d'exploration, qui permettra de limiter la profondeur de l'arbre.
  3. le coefficient de proportionnalité de la probabilité d'arrêt, qui règle en quelque sorte la vitesse d'augmentation de la profondeur d'exploration avec le temps. On le fixera à 0.9.
  4. le coefficient de proportionnalité du biais, qui permet de régler le poids accordé à l'exploration par rapport à l'exploitation.

4.3  Comparaison avec la recherche arborescente complète

Ici, nous cherchons à mettre évidence, le gain de l'UCT en terme de nombre d'états visités. Nous travaillons donc avec les paramètres fixés γ=0.9, 50 trajectoires lancées,.un coefficient pour le biais fixe, des pas de quantifications fixés à 0.02 et travaillons avec le modèle 2.2.14Modèle simplifiésubsubsection.2.2.1.

Les figures 6(a)11Subfigure 6(a)subfigure.6.1 et ?? permettent de comparer la méthode basée sur l'UCT à l'approche utilisée précédemment.
[exploration complète] [exploration par UCT avec 50 lancements de trajectoire]

Figure 6: Comparaison des deux méthodes d'exploration avec 3 actions possibles et sur une profondeur de 5.


On constate sur la figure ?? que l'on a bien quantification des états. Il faut prendre garde à la lecture de l'arbre car on n'a pas échantillonné le produit, ce qui fait que des états confondus au sens de la quantification peuvent paraître sur une ordonnée différente. Leurs abscisses sont par contre les mêmes.

On peut remarquer que les états terminaux, marqués par des points verts, les plus performants au sens de la production ne sont pas tous présents dans la figure ??. Il ne faut pas en conclure aussi vite que l'algorithme UCT fait mal son travail car ce qui nous intéresse ce n'est pas la production instantanée à un horizon temporel donné (ce qui pourrait conduire à des politiques catastrophiques), mais la somme actualisée de la production. Certes, étant donné le nombre limité de trajectoires lancées, l'algorithme UCT n'a peut-être pas trouvé la trajectoire optimale selon notre critère. Mais, on est assuré, si l'on augmente le nombre de trajectoires, de la trouver. De plus, on constate sur la figure ?? que l'arbre exploré est relativement représentatif de l'ensemble des états futurs possibles (cf. 6(a)11Subfigure 6(a)subfigure.6.1). Dès lors, on peut s'attendre à ce que le contrôle du processus, même s'il n'est pas optimal, soit performant. C'est d'autant plus vrai que la recherche arborescente par UCT va permettre d'explorer à des profondeurs bien plus grandes et donc on peut s'attendre à avoir un contrôle robuste.

Remarque: l'algorithme de recherche arborescente par UCT codé, posait un léger problème. En effet, la somme actualisée des récompenses sur un état terminal croît avec la longueur de la trajectoire. Ainsi on peut avoir une grande valeur estimée pour un état du fait qu'il fait partie d'une longue trajectoire, alors qu'une trajectoire performante mais interrompue tôt sera pénalisée par une faible valeur. Or, quand on se trouve dans la phase d'exploration, sur des états jamais visités, on a une indétermination sur l'état à choisir puisque pour toutes les actions on la même valeur (nulle) et que le biais donne le même résultat. Au départ, dans ce cas de figure, l'utilisation de la fonction max de Matlab, tranchait automatiquement le problème en choisissant la première action qui est de ne pas alimenter en substrat. Si bien que cette action était plus souvent parcourue dans l'arbre et que l'évaluation de la valeur était biaisée par cet effet. Après avoir identifié ce problème, j'y ai remédié en choisissant l'action au hasard dans les cas indéterminés.

4.4  Résultats obtenus avec le modèle 2.2.14Modèle simplifiésubsubsection.2.2.1

Maintenant que nous disposons d'une méthode de recherche arborescente un peu plus efficace, nous pouvons rechercher à reproduire l'expérience du contrôle du processus de digestion anaérobie en explorant à chaque étape sur une profondeur plus élevée.

Les paramètres que j'ai choisis sont: γ=0.9, une profondeur maximale p=10 ce qui est cohérent puisque γp=0.3487 est non négligeable devant 1, M = 4 actions, 50 itérations de l'algorithme qui suffisent à atteindre le régime permanent, un coefficient pour le biais conservé à 1, un coefficient pour l'arrêt à 0.9 et un nombre de trajectoire lancées de 60 qui permet de réaliser tous les calculs en un temps raisonnable et donc on a vérifié l'adéquation avec la profondeur par la représentation d'un arbre à partir de l'état initial. Cet état initial a été choisit de manière à ce que le régime permanent soit atteint assez rapidement.


Figure 7: Fonctionnement de l'algorithme avec l'UCT.


Le résultat du contrôle est présenté figure 712Résultats obtenus avec le modèle <A HREF="#simple">2.2.14Modèle simplifiésubsubsection.2.2.1</A>figure.7. On constate que la valeur moyenne de la production dans le régime permanent est plus élevée. On a donc su avec cette méthode d'exploration, contrôler notre système de manière plus performante. C'est étrange car la profondeur d'exploration est quasiment la même que pour l'expérience réalisée avec le modèle 2.2.14Modèle simplifiésubsubsection.2.2.1, présentée figure ??. En fait, on montre ici que pour la profondeur d'exploration considérée l'idée de choisir la première action de la trajectoire qui maximise la somme actualisée des récompenses est loin de l'action optimale.

En réalité, l'algorithme UCT codé n'est pas assez optimisé pour permettre de réellement sonder des profondeurs supérieures à celles qu'on a atteintes avec l'arborescence complète.

4.5  Résultats obtenus avec le modèle ??

Ici nous mettons en oeuvre l'algorithme pour le modèle plus pertinent mais qui fait intervenir une dimension d'espace d'état deux fois plus grande.

Bien sûr, étant connues les conclusions précédentes, nous ne nous attendons pas non plus à avoir des résultats probants. Notre prétention se limitera donc ici à présenter la mise en oeuvre.

Nous allons comparer une exploration d'arbre complète avec une exploration par UCT. Les figures 813Résultats obtenus avec le modèle <A HREF="#elabore">??</A>figure.8 présentent l'exploration complète avec 3 actions possibles et une profondeur de 5. Le nombre de noeuds de l'arbre est 364 et le temps mis à le calculer 0.090192 secondes. Les figures ?? présentent quant à elles l'exploration par UCT avec encore 3 actions possibles et une profondeur maximale de 5. Cette exploration produit 50 noeuds et prend 0.208661 secondes de calcul.


Figure 8: Recherche arborescente complète jusqu'à une profondeur de 5.




Figure 9: Recherche arborescente par UCT avec 30 trajectoires lancées et avec une profondeur maximale de 5.


On a bien vu avec le temps de calcul que l'exploration par UCT n'était pas plus performante (en termes de rapidité) pour cette profondeur d'exploration. Cependant, on constate qu'elle parvient à visiter les états les plus intéressants si l'on accorde un nombre de trajectoires suffisant.

Pour finir, nous présentons avec les figures 10(a)14Subfigure 10(a)subfigure.10.1 et ?? et une comparaison d'un contrôle simulé avec l'exploration complète sur une profondeur de 8 et avec 2 actions possibles, d'avec une exploration par UCT pour laquelle on lance 60 trajectoires et dont on vérifie bien que cela permet d'atteindre la profondeur 8.
[exploration complète] [exploration UCT]

Figure 10: Comparaison du contrôle avec une exploration jusqu'à une profondeur de 8.


Là encore, le résultat est étonnement meilleur en utilisant l'UCT.

5  Conclusion

Nous avons montré que dans le cas d'une modélisation simple du problème, un contrôle classique pouvait être largement plus aisé et performant que les méthodes par recherche arborescente. Dans un cas plus complexe, nous avons vu que la méthode de recherche arborescente était prometteuse, du fait de la faible complexité de l'arbre engendré. Cependant, le code produit au cours de ce projet ne s'est pas avéré assez performant pour pouvoir présenter des résultats probants.

Pour améliorer ce qui a été fait, on pourrait envisager de transcrire le code dans un langage compilé (par exemple C++) pour gagner en rapidité. D'autre part, au cours de l'algorithme, on pourrait réutiliser l'arbre produit à un instant pour accélérer le calcul de l'arbre suivant.

References

[1]
O. Bernard, Z. Hadj-Sadok, D. Dochain, A. Genovesi, and J.-P. Steyer. Dynamical model development and parameter identification for an anaerobic wastewater treatment process. vol. 75:424 – 438, 2001.

[2]
P. Fischer, P. Auer, and N. Cesa-Bianchi. Finite-time analysis of the multiarmed bandit problem. Machine Learning Journal, (47(2-3)):235–256, 2002.

[3]
Cs. Szepesvari and L. Kocsis. Bandit based monte-carlo planning. pages 282–293, 2006.

A  Annexes: code développé

Tout le code suivant est disponible sur ma page web: http://eleves.dptmaths.ens-cachan.fr/mguerquin/devoirs/Munos/projet.
1 % Ce script est le script principal du projet. Avec lui, on peut simuler le 2 % contrle du racteur. Les paramtres rgler sont: 3 % 1) l'tat initial via la dclaration du noeud, qui dtermine quel modle 4 % du processus utiliser par la suite. 5 % 2) le nombre d'actions M. 6 % 3) la profondeur d'exploration. 7 % 4) le nombre d'itrations de l'algorithme. 8 % 5) le coefficient d'actualisation. 9 % 6) l'algorithme utiliser pour calculer l'arbre. 10 % 11 % A la fin, on affiche la trajectoire suivie au cours du temps dans 12 % l'espace d'tat. 13 14 15 M = 2; % nombre d'actions 16 p = 8; % profondeur 17 nbr = 100; % itrations 18 gamma = 0.9; % coefficient d'actualisation 19 20 % structure de chaque noeud + intialisation 21 %noeud = declaration_noeud([1 1 0],M); %modele1 22 %noeud = declaration_noeud([0.4 6 0],M); %modele1[1 1 0] 23 noeud = declaration_noeud([1 1 1 1 0],M); %modele2 24 25 % Modle utiliser 26 modele = @modele1; 27 if length([fieldnames(noeud.etat)])>3 28 modele = @modele2; 29 end 30 31 % Dbut de la simulation 32 data = []; 33 action = []; 34 for i = 1:nbr 35 % sauvegarde de la trajectoire 36 data = [data noeud.etat]; 37 % Calcul de l'arbre 38 arbre = exploration(noeud,[1:M],gamma,p,M); 39 %arbre = UCT(noeud,[1:M],gamma,p,M); 40 41 % Choix de l'action optimale 42 [somme_act sequence] = best_somme_act(arbre); 43 44 45 % Action constante 46 %sequence(1) = 2; 47 %arbre(1) = noeud; 48 49 50 51 action = [action sequence(1)]; 52 % calcul de l'tat suivant 53 noeud.etat = modele(arbre(1).etat,sequence(1),M); 54 end 55 56 %%%%%%%%%%%%%%%%%%%%%%%%%%% 57 % Affichage des rsultats 58 if length([fieldnames(noeud.etat)])<4 59 figure(1); 60 plot3([data.population1],[data.substrat1],[data.produit],'b*-'); 61 xlabel('population');ylabel('substrat');zlabel('produit'); 62 hold on; 63 h = plot3([data(1).population1],[data(1).substrat1],[data(1).produit],'r.'); 64 set(h, 'MarkerSize', 25); 65 h = plot3([data(end).population1],[data(end).substrat1],[data(end).produit],'g.'); 66 set(h, 'MarkerSize', 25); 67 hold off; 68 legend('trajectoire','point de dpart','point d''arrive',0); 69 end 70 71 figure(2); 72 x = 1:length([data.population1]); 73 plot(x,[data.population1],'r-*');hold on; 74 plot(x,[data.substrat1],'g-+'); 75 plot(x,10*[data.produit],'b-<'); 76 plot(x,[action-1]/(M-1),'k-*'); 77 if length([fieldnames(noeud.etat)])>3 78 plot(x,[data.population2],'c->'); 79 plot(x,[data.substrat2],'m-o'); 80 legend('population1','substrat1','10*production','action','population2','substrat2',0); 81 else 82 legend('population1','substrat1','10*production','action',0); 83 end 84 hold off 85 xlabel('temps');

1 function arborescence(etatinit,M,p) 2 % Cette fonction sert lancer l'affichage d'un arbre partir d'un tat 3 % initial. 4 % M: nombre d'actions 5 % p: profondeur de l'arbre 6 % etatinit: tat initial: [sub1 pop1 prod] ou [sub1 pop1 sub2 pop2 prod] 7 8 % arborescence([0.32,0.4,0.02],4,3) 9 % arborescence([0.36,0.35,0.15,0.16,2e-3],4,3) 10 11 % structure de chaque noeud 12 noeud = declaration_noeud(etatinit,M); 13 14 % Calcul de l'arbre 15 tic;arbre = exploration(noeud,[1:M],0.9,p,M);toc 16 %tic;arbre = UCT(noeud,[1:M],1.2,p,M);toc 17 18 length(arbre) 19 20 % Affichage 21 plot_arbre(arbre);

1 function arbre = exploration(noeud,actions,gamma,p,M) 2 % Explore, de manire nave, l'arbre des tats possibles partir de 3 % 'noeud' et en suivant les M actions, en s'arrtant une profondeur p. 4 % La somme des rcompenses est actualise avec le coefficient gamma. 5 6 arbre = [noeud]; 7 % Quel modle utiliser? 8 modele = @modele1; 9 if length([fieldnames(noeud.etat)])>3 10 modele = @modele2; 11 end 12 13 for i = 1:p %pour chaque profondeur 14 for x = search_depth(arbre,i-1) %pour chaque noeud la profondeur inf 15 for a = actions % pour chaque action 16 % calcul du nouvel tat 17 noeud.etat = modele(arbre(x).etat,a,M); 18 % mise jour des donnes 19 noeud.profondeur = i; 20 noeud.parent = x; 21 noeud.suivant = []; 22 arbre(x).somme_act(a) = arbre(x).somme_act(a)+gamma^(i-1)*noeud.etat.produit; 23 noeud.somme_act(:) = arbre(x).somme_act(a); 24 arbre(x).suivant(a) = length(arbre)+1; 25 arbre = [arbre noeud]; 26 end 27 end 28 end

1 function arbre = UCT(noeud,actions,gamma,p,M) 2 % Explore, en utilisant l'algorithme UCT, l'arbre des tats possibles 3 % partir de 'noeud' et en suivant les M actions. Les rcompenses sont 4 % actualises avec le coefficient gamma. L'arbre va jusqu' la profondeur 5 % maximale p. 6 % 7 % Le nombre de trajectoires lances (timeout), 8 % le facteur de proportionnalit de la probabilit d'arrt, 9 % et le facteur de proportionnalit du biais sont ajuster ci dessous. 10 11 12 timeout = 60; 13 coeff_stop = 0.7; 14 coeff_biais = 0.1; 15 16 arbre = [noeud]; 17 % Quel modle utiliser? 18 modele = @modele1; 19 if length([fieldnames(noeud.etat)])>3 20 modele = @modele2; 21 end 22 23 for t = 1:timeout 24 %lancement d'une trajectoire partir de l'tat initial 25 arret = 0; 26 x = 1; % indice du noeud courant 27 arbre(x).compteur_etat = arbre(x).compteur_etat+1; 28 while ~arret % tant qu'il n'y a pas de condition d'arret 29 % choix de l'action 30 %arbre(x).somme_act(:)+bias_term(arbre(x)) 31 valeur = estimate_value(arbre,x); 32 biais = bias_term(arbre(x),coeff_biais); 33 [tmp a] = max(valeur(:)+biais(:)); 34 test = find((valeur(:)+biais(:))==tmp); 35 if length(test)>1 36 a = test(round((length(test)-1)*rand(1))+1); 37 end 38 % calcul du nouvel tat chantillonn si besoin 39 if (arbre(x).suivant(a)==0) %si suivant pas encore visit 40 % on le calcule 41 noeud.etat = echantillonnage(modele(arbre(x).etat,a,M)); 42 noeud.profondeur = arbre(x).profondeur+1; 43 noeud.parent = x; 44 % on vrifie si cet tat existe dj dans l'arbre 45 y = search_depth(arbre,noeud.profondeur); 46 if (length(y)~=0) 47 visited = find(egalite_etat([arbre(y).etat],noeud.etat)); 48 visited = y(visited); 49 else 50 visited =[]; 51 end 52 if (length(visited)~=0) % Si l'tat est dj dans l'arbre 53 %%%%%%%%% 54 % CAS 1 55 %on actualise les somme_acts et le parent du noeud 56 arbre(x).somme_act(a) = arbre(x).somme_act(a)+gamma^(arbre(x).profondeur-1)*arbre(visited).etat.produit; 57 % Changement du parent de l'tat suivant? 58 ind = arbre(visited).parent; 59 % C'est une question de somme actualise ! 60 if arbre(ind).somme_act(find([arbre(ind).suivant]==visited))<arbre(x).somme_act(a) 61 arbre(visited).parent = x; %changement de parent 62 arbre(visited).somme_act(:)=arbre(x).somme_act(a); %et de somme_act 63 end 64 arbre(x).suivant(a) = visited; 65 else %si l'tat suivant n'est pas dans l'arbre 66 %%%%%%%%%% 67 % CAS 3 68 arbre = [arbre noeud];% on l'y ajoute 69 arbre(x).somme_act(a) = arbre(x).somme_act(a)+gamma^(arbre(x).profondeur-1)*noeud.etat.produit; 70 arbre(x).suivant(a) = length(arbre); 71 arbre(end).somme_act(:) = arbre(x).somme_act(a); 72 end 73 else % si l'tat suivant est dj visit 74 %%%%%%%%%% 75 % CAS 2 76 end 77 %on actualise les compteurs et le parent du noeud 78 arbre(arbre(x).suivant(a)).compteur_etat=arbre(arbre(x).suivant(a)).compteur_etat+1; 79 arbre(x).compteur_etatactions(a)=arbre(x).compteur_etatactions(a)+1; 80 81 % Arrte-t-on l la trajectoire? 82 if (arbre(x).profondeur == p-1)|(rand(1)<coeff_stop/(1+0.5*arbre(x).compteur_etat)) 83 arret==1; 84 break; 85 end 86 x = arbre(x).suivant(a); 87 end 88 end 89 90 function bool = egalite_etat(etat1,etat2) 91 % teste l'galit de deux tats 92 93 if nargin < 2 94 bool = 0; 95 else 96 bool = ([etat1(:).population1]==etat2.population1).*([etat1(:).substrat1]==etat2.substrat1); 97 if length([fieldnames(etat1(1))])>3 98 bool = bool.*([etat1(:).population2]==etat2.population2).*([etat1(:).substrat2]==etat2.substrat2); 99 end 100 end 101 102 function C = bias_term(noeud,coeff_biais) 103 % Calcule le biais a partir des nombres de visites de l'tat et de ses 104 % tat-action 105 106 C = 2*coeff_biais*sqrt(log(noeud.compteur_etatactions(:)+1.01)./(noeud.compteur_etatactions(:)+0.01));

1 function etat_suiv = modele1(etat,action,M) 2 % Calcule l'tat suivant en utilisant le modle de connaissances et 3 % l'action choisie. (M sert normaliser l'action entre 0 et 1) 4 5 % Dfinition des paramtres 6 alpha = 0.1; 7 k1 = 1; 8 k2 = 2; 9 k3 = 2; 10 c1 = 1; 11 c2 = 3; 12 c3 = 8; 13 Sin = 0.5; 14 etat_suiv = etat; 15 action = (action-1)/(M-1); 16 tau = 1; 17 18 etat_suiv.substrat1 = etat.substrat1 +tau*(action*(Sin - etat.substrat1)-k1*etat.population1*c1*etat.substrat1/(c2*etat.substrat1^2+c3)); 19 if etat_suiv.substrat1<0 20 etat_suiv.substrat1=0; 21 end 22 etat_suiv.population1 = etat.population1 +tau*(-alpha*action*etat.population1+k2*etat.population1*c1*etat.substrat1/(c2*etat.substrat1^2+c3)); 23 if etat_suiv.population1<0 24 etat_suiv.population1=0; 25 end 26 etat_suiv.produit = k3*etat.population1*c1*etat.substrat1/(c2*etat.substrat1^2+c3);

1 function etat_suiv = modele2(etat,action,M) 2 % Calcule l'tat suivant en utilisant le modle de connaissances 2 et 3 % l'action choisie. (M sert normaliser l'action entre 0 et 1) 4 5 % Dfinition des paramtres 6 alpha = 0.2; 7 k1 = 1; 8 k2 = 2; 9 k3 = 2; 10 k4 = 2; 11 c1 = 1; 12 c2 = 8; 13 c3 = 1; 14 c4 = 3; 15 c5 = 8; 16 S1in = 0.5; 17 S2in = 0.1; 18 etat_suiv = etat; 19 action = (action-1)/(M-1); 20 tau = 1; 21 22 etat_suiv.substrat1 = etat.substrat1 +tau*(action*(S1in - etat.substrat1)-k1*etat.population1*c1*etat.substrat1/(c2+etat.substrat1)); 23 if etat_suiv.substrat1<0 24 etat_suiv.substrat1=0; 25 end 26 etat_suiv.population1 = etat.population1 +tau*(-alpha*action*etat.population1+etat.population1*c1*etat.substrat1/(c2+etat.substrat1)); 27 if etat_suiv.population1<0 28 etat_suiv.population1=0; 29 end 30 etat_suiv.substrat2 = etat.substrat2 +tau*(action*(S2in - etat.substrat2)+k2*etat.population1*c1*etat.substrat1/(c2+etat.substrat1)-k3*etat.population2*c3*etat.substrat2/(c4*etat.substrat2^2+c5)); 31 if etat_suiv.substrat2<0 32 etat_suiv.substrat2=0; 33 end 34 etat_suiv.population2 = etat.population2 +tau*(-alpha*action*etat.population2+etat.population2*c3*etat.substrat2/(c4*etat.substrat2^2+c5)); 35 if etat_suiv.population2<0 36 etat_suiv.population2=0; 37 end 38 etat_suiv.produit = k4*etat.population2*c3*etat.substrat2/(c4*etat.substrat2^2+c5);

1 function noeud = declaration_noeud(etatinit,M) 2 % Fonction permettant de dclarer un noeud pour l'arbre en l'initialisant 3 % l'tat dfinit par 'etatinit'. 4 5 if length(etatinit)==3 6 noeud.etat = struct('substrat1',etatinit(1),'population1',etatinit(2),'produit',etatinit(3)); 7 else 8 noeud.etat = struct('substrat1',etatinit(1),'population1',etatinit(2),'substrat2',etatinit(3),'population2',etatinit(4),'produit',etatinit(5)); 9 end 10 noeud.profondeur = 0; 11 noeud.parent = [0]; %indice du parent 12 noeud.suivant = [zeros(1,M)]; %indice des enfants selon l'action 13 noeud.somme_act = [zeros(1,M)]; %somme actualise des rcompenses associe chaque action 14 noeud.compteur_etat = 0; %initialisation du compteur de l'tat 15 noeud.compteur_etatactions = [zeros(1,M)]; %initialisation des compteurs

1 function etat = echantillonnage(etat) 2 % fonction ralisant l'chantillonnage des tats avec les pas ci-dessous 3 4 pas_pop = 0.02; 5 pas_sub = 0.02; 6 7 etat.substrat1 = pas_sub*round(etat.substrat1/pas_sub); 8 etat.population1 = pas_pop*round(etat.population1/pas_pop); 9 if length([fieldnames(etat)])>3 10 etat.substrat2 = pas_sub*round(etat.substrat2/pas_sub); 11 etat.population2 = pas_pop*round(etat.population2/pas_pop); 12 end

1 function plot_arbre(arbre) 2 % Fonction qui permet de reprsenter un arbre dans l'espace d'tat 3 4 close all; 5 6 racine = arbre(1); 7 M = size(racine.somme_act); %nombre d'actions 8 p = max([arbre.profondeur]); %profondeur de l'arbre 9 etats = [fieldnames(racine.etat)]; 10 n = length(etats); 11 color = ['b' 'g' 'r' 'c' 'm' 'y' 'k']; 12 13 for i = 1:p 14 ind = search_depth(arbre,i); 15 for j = ind 16 col = color(mod(getfield(arbre(j),'profondeur')-1,length(color))+1); 17 for fig = 1:n-1 18 figure(fig); 19 hold on; 20 nom = char(etats(fig)); 21 prod = char(etats(end)); 22 plot([getfield(arbre(arbre(j).parent).etat,nom) getfield(arbre(j).etat,nom)],[getfield(arbre(arbre(j).parent).etat,prod) getfield(arbre(j).etat,prod)],[col '-*']); 23 hold off; 24 end 25 end 26 end 27 for fig = 1:n-1 28 figure(fig) 29 hold on; 30 h = plot([getfield(arbre(1).etat,char(etats(fig)))],[getfield(arbre(1).etat,char(etats(end)))],'r.'); 31 set(h, 'MarkerSize', 25); 32 ind = search_leaves(arbre); 33 for i = ind 34 h = plot([getfield(arbre(i).etat,char(etats(fig)))],[getfield(arbre(i).etat,char(etats(end)))],'g.'); 35 set(h, 'MarkerSize', 25); 36 end 37 hold off; 38 xlabel(char(etats(fig))); 39 ylabel(char(etats(end))); 40 end

1 function [somme_act sequence] = best_somme_act(arbre) 2 % Fonction qui cherche le noeud parmis les feuilles de l'arbre qui possde 3 % la plus grande somme actualis et retourne la suence d'actions qui 4 % mnent ce noeud ainsi que sa somme actualise. 5 6 ind = search_leaves(arbre); 7 a=[arbre(ind).somme_act]; 8 b=a(1:length(arbre(1).somme_act):end); 9 [somme_act x] = max(b); 10 sequence = []; 11 ind = ind(x); 12 for i=1:arbre(ind).profondeur 13 ind2 = arbre(ind).parent; 14 sequence = [find([arbre(ind2).suivant]==ind) sequence]; 15 ind = ind2; 16 end

1 function valeur = estimate_value(arbre,ind_noeud) 2 % Fonction rcursive qui estime la valeur d'un noeud dans l'arbre en 3 % cherchant parmi les feuilles qui en sont issues la somme actualise 4 % maximale. 5 6 nbr = length(ind_noeud); 7 if nbr>1 8 for i=1:nbr 9 valeur(i) = estimate_value(arbre,ind_noeud(i)); 10 end 11 end 12 13 s = sum([arbre(ind_noeud).suivant]); 14 15 if (s==0) 16 valeur = arbre(ind_noeud).somme_act; 17 else 18 i = find(arbre(ind_noeud).suivant(:)~=0); 19 seq = []; 20 for j = 1:length(i) 21 seq = [seq estimate_value(arbre,arbre(ind_noeud).suivant(i(j)))]; 22 end 23 valeur = max(seq); 24 end

1 function ind = search_depth(arbre,p) 2 % retourne les indices des noeuds de l'arbre qui sont la profondeur p 3 4 ind = find([arbre.profondeur]==p);

1 function ind = search_leaves(arbre) 2 % Fonction qui retourne les indices des feuilles de l'arbre 3 4 s = zeros(1,length(arbre)); 5 for i = 1:length(arbre) 6 s(i) = sum(arbre(i).suivant(:)); 7 end 8 ind = find(s

1
mailto:matthieu.guerquin-kern AT ens-cachan DOT fr
2
sources: http://www.ercim.org/telemac/pub/factsheet-F.pdf
This document was translated from LATEX by HEVEA.